General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



FACILITY FORM 602 



THORIZATION 


Control No . 
Date _____ _ 


MSC INTERNAL NOTE 


EFFICIENCY OF GENERALIZED 
MATRIX INVERSION METHODS 


Fred C. Delaney 


C ary G. Gaffney 


Fred M 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATE) N 
MANNED SPACECRAFT CENTER 
HOUSTON, TEXAS 

September 1, 1966 


(ACCESSION NUMBER) 


(PAGES) 


(CATEGORY) 


(NASA CR OR TMX OR AD NUMBER) 


MSC INTERNAL NOTE NO. 66-ED- \l 


Prepared by 


Prepared by 


Prenared by 


Approved by 


Approved by 


NATIONAL 


EFFICIENCY OF GENERALIZED 
MATRIX INVERSION METHODS 


Fred C. Df mey 

Lockheed Electronics Company 



Aerospace Writing and Analysis Program 

T n. JLu l . 

Fred M. Speecf^ / 

Theory and Analysis Office 





L. Davis, Jr. 1 / 

Theory and Analysis Office 

L 0^t jL — - 

Eugerje h . Brock 


Eugei 

Chi ey, Computation and 
Analysis Division 


AERONAUTICS AND SPACE ADMINISTRATION 
MANNED SPACECRAFT CENTER 
HOUSTON, TEXAS 

September 1, 1966 


TABLE OP CONTENTS 


Page 

SUMMARY , 1 

INTRODUCTION 1 

. - 

SYMBOLS • ' ■ 2 

METHODS STUDIED • 

* 

• Accuracy Determinations - Methods and Results 6 

Results of Speed Determinations 9 

CONCLUSION 11 

i 

APPENDIXES 

A. Theoretical Background ' A~1 

B. Analysis of Variance ' ‘ B-l 


REFERENCES 


LIST OP ILLUSTRATIONS 


Figure 

1 Graph of Time Vs . Order for GINV2 

* , 

LIST OP TABLES 

» » 

Table 

I Norm Values for GINV2 and APLUS for 

Selected Randomly Generated Matrices 

II Computer Times Required to Calculate the 

Generalized Inverse of a 10x10 Non~ 

• singular • Matrix 

III Computer Times Required by GINV2 and 

APLUS for Selected Randomly Generated 
Matrices 

IV Computer Times Required by GINV2 and 

APLUS for Selected Randomly Generated 
Matrices 


Page 

12 

* 

Fage 

13 

21 

21 

22 


EFFICIENCY OF .GENERALIZED 
MATRIX INVERSION METHODS 

. •• 

By Fred C. ( Delaney, Gary G, Gaffney, and Fred M. Speed 

. SUMMARY 

t 

Generalized matrix inversion has become an increasingly 
important concept in matrix theory as well as a useful tool 
in engineering, statistics, control theory, and space mission 
design. For this reason, the need arises for an efficient 
(i.e., fast and accurate) method of computing the general- 
ized inverse of an arbitrary n x m matrix. The purpose of 
this paper is to present the results of computer tests used 
to compare the relative efficiency of several computer pro- 
grams designed to calculate the generalized inverse of 

an arbitrary real matrix. 

* * 

•> 

INTRODUCTION 

The concept of matrix inversion was first generalized 
by E. H. Moore (ref. 1) in 1920. In the 1950’s 
R. Penrose (ref. 2) and A. Bjerhammer (ref. 3), working 
independently, formulated equivalent definitions of the 
generalized inverse of an arbitrary complex matrix. 

The most common definition, given by Penrose (ref. 2) is 
a consequence of the following, theorem. .* 



THEOREM 1: For any real n x m matrix A, there is a unique 

real m x n matrix A* (the generalized inverse of A) such 
that : 

(1) AA^A = A 

(2) A + AA + = A + ' 

» 

(3) ’ •(AA + ) T = AA* 

(4) (A + A) T = A + A 


The use of the generalized inverse in engineering 
problems, statistics, and control theory gave rise, 
naturally enough, to the development of several different 
computational methods. Some of these methods were developed 
and programed by researchers at NASA-MSC to solve problems 
requiring generalized matrix inversion. 

* * 

This paper presents the results of an examination of 

the various programs for overall efficiency, comparing 
them in terms of ’ accuracy and speed. 

SYMBOLS 

Capital letters denote matrices with real entries. 

a n is a row vector. . . • 

is a column vector. 

n 

T 

A denotes the transpose of the matrix A. 


•fk-ft*! 


*•*- 






} „ •*»*»** *»«**;. 


A + denotes the generalized inverse of the matrix A. 

th th 

(A)ij is the entry in the i row and j column of the 
matrix A. 

• » 

» * 

% * 

|| || is a matrix norm. 

• » 

Trace (A) is the trace of the matrix A. ( 

* * 
t 

I is the identity matrix. 

1^ is the k x k identity matrix. 

(A;B) is a matrix partitioned into the matrix A and the 

matrix B. 

, * 


t 










METHODS STUDIED 


Although there are numerous mathematical methods for 
calculating the generalized Inverse of a matrix, the purpose 
of this study was to determine the most efficient (that Is, 
fastest and most accurate) computer method for calculating 
this Inverse for arbitrary real matrices, A preliminary 
survey of .the existing algorithms for generalized matrix 
inversion showed that some of them were not readily adapt- 
able to computer programing or were more suitable only to 
theoretical investigations and required no further cons id.-' 
eration. One of the Penrose methods (ref. k) was discarded 
because it first requires a type of matrix partitioning 
that is time consuming on the computer. The Ben-Israel and 
Wersan method (ref. 5) was eliminated because it depends on 
the exact determination of rank, which depends on round-off 
and approximation errors. The Householder method (ref. 6) 
was rejected because it depends on predetermining the rank 
of the matrix. Since the Ben-Israel and Charnes method 
(ref. 7) uses the Lagrange-Sylvester interpolation poly- 
nomial, which is sensitive to error in the computer, it, too 
was discarded. Finally, the Decell method based on the 
Cayley-Hamilton theorem (ref. 8), which requires the cal- 
culation of powers of a matrix, was eliminated because of 
the error which such a calculation causes . 

The following five methods, having satisfied this pre- 
liminary requirement of computer adaptability, were then 
examined because they showed promise of being efficient 


4 




generalized Inverse programs: 

7 r -9 * • * musmttx* # •«•** ■ - * #-*-#«»-•*«: - -rt -••» . 2i R , . at , 

(1) the computer program PEN2, based on another of the 
Penrose (ref, 2) methods , 

(2) the program SEQ1NV, based on a method by 
H. P. Decell, Jr., (ref, 9) of NASA-MSC, 

(3) JPLUS, a method developed and programed by two of 
the authors, F. M. Speed of NASA-MSC, and 

F« C, Delaney of DEC, 

(4) A PLUS, taken from an iterative method devised by; 
H. P„ Decell, Jr,, and $. W, Kahng (ref, 10), 

(5) GINV2 , an algorithm developed by B. Rust, 

» 

W. H, Burras, and C* Schneeberger (ref. 11). 

*v 

(For the mathematics underlying these methods see 
Appendix A.) 

These methods were programed in FORTRAN IV, if they had 
not already been, and were then tested very extensively on 

the UNIVAC 1108 to determine efficiency. 

\ 

Since generalized matrix inversion is applicable to 
arbitrary matrices, some preliminary mention should be made 
of the variety of matrices used in testing the programs . 

The generalized inverses of singular and nonsingular square 
matrices and of nonsquare matrices of full and less than 
full rank were computed by each of the five methods. These 
matrices were, for the most part, randomly generated and 
differed in size from order 2 x 2 to order 45 x 40. As the 
results will indicate, the type and size of the matrix, 
whether due to round-off error in the computer, or to the 


increased computer time required by larger matrices, or to 
peculiarities in the method of generalised inversion, often 
had significant effects on the speed and accuracy of the 
program. 

* 

* 

Accuracy Determinations - Methods and Results 

Before any test results can be presented, a description 
of the methods for determining and comparing the accuracies 
of the above five programs is necessary. 


The four identities of THEOREM 1, which define the gen- 
eralized inverse, suggest a means for testing the accuracy 
of a program designed to calculate it. In the case of real 
matrices, norms based on these identities can be defined in 
the following way : 


Let A be an n x m matrix with. real entries, and let A . 
denote the generalized inverse ’as calculated by computer. 

A 

Then A is an m x n matrix also with real entries. 


Define : 


NORM 1 = j | AAA -A|| 


NORM 2 = | | AAA - A| | 



nm 



« 



These norms provide a satisfactory test for accuracy 
since, first, each is, in fact, the root mean square of an 
element in the difference matrix for that norm and, second, 

* j. 

each norm is equal to zero if and only if A is equal to A , 
the true generalized inverse of A. Note, however, that these 
norms will not, as a rule, be zero due to round-off error 

in the computer. 

. / •• 


in order for a generalized inversion program to be con- 
sidered dependable and to have wide application, it must be' 
capable of computing the generalized inverse of all types 

of matrices with a consistent and predictable accuracy. 

* 

The results of the tests in this study showed that not all 
of the five programs mentioned above could meet this demand. 

Two programs, APLUS and GINV2, did, however, perform 
with more than satisfactory accuracy for all of the matrix 
types used in testing. The norms evaluated using the v 

generalized inverses computed by these two programs ranged 
nearly always between 1 x 10” 4 and 1 x 10~ 12 , averaging 
between 1 X 10 and 1 x 10 . (For a more detailed com- 

parison see Table 1.) 


In order to determine if there was a significant 
difference in accuracy between APLUS and GINV2, an analysis 
of variance was performed (see Appendix B for the details 

of the analysis). The results of the analysis showed that, 

» 

for full rank matrices, it is not possible to reject the 
hypothesis that the accuracies of the two methods are equal. 
However, for matrices of nonfull rank (excluding n x m 
matrices of rank 1), the hypothesis that the accuracies of 
the two methods were equal was rejected in favor of the 
hypothesis that GINV2 was more accurate than APLUS . 

The remaining three methods, PEN2, SEQINV, and JPLUS, 

could not consistently meet the demands on accuracy, and 

» • 

therefore will certainly have restrictions~~in varying 
degrees — on their application. ' 

PEN2 gives acceptable norms for some small matrices 
and for matrices of very low rank but gives very poor results 
for all other types of matrices . This program is not very 
dependable and should find little, if any, application. 

SEQINV yields good results for nonsingular matrices 
and matrices with full rank or low rank, but as the order 
of the matrix increases past 30 x 30 it begins to fail 
noticeably for singular matrices and matrices with less than 
full rank. Even when SEQINV performs well, its accuracy does 
not exceed — and usually lags behind-.--that of APLUS and GINV2 . 
For completeness, it should be noted here that SEQINV con- 
tains a zero-test whose epsilon value, when increased slightly, 
causes significant improvement in some norms which had 
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previously indicated that the program had failed, This 

epsilon value was not experimented with in detail since 

varying it caused no significant change in those norms for 

which SEQINV , formerly gave, good results. 

% * 

Of the above three methods, JPLUS is, by far, the most 
consistently' accurate and dependable. It yields poor norms 
in only a small number of cases j namely where the matrices 
were singular and of large rank and order, hut despite its 
rather satisfactory performance, its accuracy is not as 
great as that of APLUS or GINV2, which limits its application. 
JPLUS, like SEQINV, also has a zero-test, whose epsilon 
value, when varied, .causes changes in the norms with results 
very similar to those observed for SEQINV. 

■v 

Results of Speed Determination 

• In order to obtain a sample of the relative speed of 
each program, a test; block consisting of one hundred 10 x 10 
nonsingular matrices was generated randomly. The computer 
time required by each program (except PEN2) for calculating 
£he generalized inverse of each matrix in the block was then 
determined for comparison purposes. (See Table II.) Because 
the levels of accuracy for PEN2, SEQINV, and JPLUS were not 
entirely satisfactory, no further time tests were made on 
these programs. 

Since APLUS and GINV2 were the only programs which met 
the demands on accuracy, and since they were found to'have 
nearly equal accuracy, time was the deciding factor in 


determining which was the most efficient computer program 
examined. For this reason, time samplings were run for 
blocks of 20 x 20, 30 x 30, and *10 x *10 nonsingular 
matrices of the type described above. (To conserve computer 
time, the number of matrices per block was reduced as the 
matrix order increased.) The results of these samplings 
showed that GINV2 is considerably faster than APLUS. (See 
Table III..) 

The GINV2 program must do an additional set of operations 

for any dependent column in a matrix whose generalized 

inverse is to be computed. For this reason, the times 

required by APLUS and GINV2 were also compared for singular 

♦ 

matrices. Various blocks, each of 15 matrices of the same 
rank and order, were again generated randonly. The matrix 
types tested were order 10 x 10 matrices of ranks 1 through 
10 and order 20 x 20 matrices of ranks 1 through 20. (See 
Table IV.) 

t 

It was observed that for matrices of rank 1, both of 
order 10 x 10 and of order 20 x-20, APLUS is slightly faster 
than GINV2. It seerns reasonable to conclude that for the 

v 

rank 1 case APLUS makes a sufficiently accurate initial guess 
at the generalized inverse and the iteration process stops 
immediately. However, for matrices of rank 2, APLUS 
becomes considerably slower while GINV2 becomes slightly 
-faster so that the times for the two methods compare much 
as they did in the nonsingular case. As the rank increases 
to full rank, APLUS becomes generally slower while GINV2 
becomes increasingly faster. For a fixed matrix size, GINV2 



is at its fastest when the rank is maximum. It was observed 
that GINV2 averages 10 to 15 times faster than APLUS . 

CONCLUSION 

% • % 

The results of this study clearly indicate that GINV2 

is the most efficient program (among the computer subroutines 

studied) for calculating the geheralized inverse of a matrix. 

(See Figure 1.) Both APLUS and GINV2 are dependable methods 

in terms of accuracy, but GINV2 is considerably faster. 

APLUS is more efficient than GINV2 in only one special 

case — matrices (other than m x 1 matrices) of rank 1; and 

* 

In this case the following simple formula exists for com- 
puting the generalized Inverse: 

% 

A+ 1 «T 

A . ** m A • 

trace(A A) 

Further information as well as copies of the computer 
programs can be obtained from: 

i 

. F. M. Speed 

, 'Theory, and Analysis Office 

National Aeronautics and Space Administration 
Manned Spacecraft Center 
Houston, Texas 
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FIGURE 1 

GRAPH OF TIME VS. ORDER FOR GINV2 





Table I 

Norm Values For GINV2 and APLUS 
For Selected Randomly Generated Matrices 
(The norm values for APLUS appear first for each matrix type) 


ORDER RANK 


x 5 


x 5 


H x 5 










NORM 1 

NORM 2 . 

NORM 3 

NORM 1 

2.9 x 10" 8 

1.1 x 10“ 10 

0 

0 

1 . 1 ) x 10“ 7 

5.7 x 10" 30 

2.6 x 10“ 9 

2.6 x 10" 8 

1.6 x lCf 6 

8.7 x 10" 9 

< 1.0 x 10 -9 

< i.O x 10" 9 

i 

o 
1 — 1 

X 

oo 

• 

on 

1.1 x 10" 9 

0 

0 

< 1.5 x 10~ 7 

5.5 x 10 " n 

0 

2.2 x 10" 9 

2 . i ) x 10" 8 

1.2 x 10" 13 

0 

< 1.3 x 10" 8 

< l.<l x 10" 6 

1.7 x 10 “'’ 

< 1.3 X 10" 9 

1.0 x 10 * 8 

1.2 x 10" 7 

2.0 x 10" 9 

1 .] x 10" 9 

9.9 x 10~ 9 

1.6 x 10" 6 

1.7 x 10" 31 

6.6 x 10" 30 

' 0 

1.5 x 10" 6 

3.3 x 10" 11 

1.3 x 10" 9 

6.9 x 10" 8 

3.0 x 10” 7 

1.1 x 10" 8 

1.3 x 10" 8 

1.2 x 10* 8 

l .<) x 10" 6 

< 1.2 x 10" 10 

6.6 x 10 -9 

1.7 x 10" 7 

2.1 x 10" 7 

8 . <1 x 10" 10 

7.3 x 10" 9 

1.1 x 10" 8 

2.2 x 10“ 6 

6.0 x 10" 9 

2 . <) x 10" 7 

1. 3 x 10" 7 

2.2 x 10" 6 

.... hm 

2.9 x 10 -13 

7.1 x 10" 10 

9.2 x 10" 10 

< 1.6 x 10" 6 

8.1 x 10" 11 

9.2 x 10" 10 

< 4.8 X 10" 7 
















































Table I (Continued) 


ORDER RANK 


10 x 8 2 


10 x 8 5 


10 x 8 8 


10 x 10 



8 x 10 2 


8 x 10 5 


8 x 10 8 


10 x 10 2 


NORM 1 

NORM 2 

NORM 3 

NORM 4 

3.5 

X 

10- 6 

5.4 

X 
• • 

10“ 9 

1.4 

X 

10- 8 

4.1 

X 

10- 8 

3.9 

X 

10- 6 

2.6 

X 

10- 10 

• 

• 

7.1 

X 

10~ 9 

1.2 

X 

10~ 7 

7.6 

X 

10-7 

1.0 

X 

10- 8 

1.7 

X 

:o - 8 

1.7 

X 

10- 8 

1.7 

X 

10- 6 

3.0 

X 

y-> 

O 

1 

y-> 

O 

6.7 

X 

] 0~ 9 

2.1 

X 

10~ 7 

1.0 

X 

1 C ' 6 

2.1 

X 

ID ' 8 

1.4 

X 

10- 8 

1.8 

X 

10- 8 

4.2 

X 

] 0" 7 

9-6 

X 

IQ ' 10 

1.3 

X 

10- 8 

1.5 

X 

10- 8 

2.9 

X 

10" 6 

2.7 

X 

10- 11 

5.4 

X 

] 0- 10 

4.5 

X 

V-* 

o 

1 

M 

o 

5.7 

X 

10- 6 

4.8 

X 

10 -11 

1.1 

X 

10” 9 

3.8 

X 

10- 7 

2.9 

X 

10- 6 

3.9 

X 

10~ 9 

9.5 

X 

10” 9 

1.1 

X 

10- 8 

• 

CD 

X 

10- 6 

8.4 

X 

10- 11 

4.2 

X 

] 0'* 9 

6.9 

X 

10” 7 

• 

OO 

\ 

X 

10- 6 

3.1 

X 

10- 8 

1.4 

X 

10- 8 

2 .,7 

X 

3 . 0- 1 8 

1.4 

i. 

X 

10" 6 

6.2 

X 

10- 10 

7.2 

X 

10" 9 

3.5 

X 

10" 7 

3.4 

X 

10~ 7 

1.0 

X 

10~ 9 

9.4 

X 

10" 9 

2.0 

X 

10- 8 

5.5 

X 

1 (T 7 

1.2 

X 

10~ 9 

2.8 

X 

10- 8 

6.2 

X 

10- 8 

3.1 

X 

10- 6 

2.0 

X 

10- 11 

5.4 

X 

10- 10 

4.5 

X 

io- ]0 

9.4 

X 

10~ 6 

5 . 7 

X 

lO” 11 

1.3 

X 

10” 9 

3.8 

X 

10“ 7 

5.9 

X 

10- 6 

8.3 

X 

10~ 9 

1.2 

X 

10- 8 

3.4 

X 

10- 8 1 

7.3 

X 

10- 6 

3.2 

X 

io' 10 

1.7 

X 

10” 9 

4.5 

X 

10“ 7 

























































Table I (Continued) 


ORDER RANK 


10 x 10 5 


30 x 10 


10 x 10 10 



20x20 2 



20 x 20 5 



20 x 20 10 


20 x 20 15 


20 x 20 20 



NORM 1 

NORM 2 

NORM 3 

NORM *1 

• 

oo 

X 

10- 6 

2.6 

X 

10- 8 

1.4 

X 

10- 8 

1.3 

X 

10~ 7 

3.1 

X 

10- 6 

3.4 

X 

10- 10 

5.9 

X 

10" 9 

2.9 

X 

10~ 7 

1.0 

X 

10- 6 

1.6 

X 

10- 8 ' 

4.1 

X 

10- 8 

1.7 

X 

10~ 7 

9.0 

X 

10~ 7 

1.1 

X 

10- 8 

5.6 

X 

10~ 8 

SK 6 

X 

10- 8 

e.n 

X 

10~ 5 

1.8 

X 

10” 11 

9.0 

X 

ht ™ 

4.5 

X 

10- 10 

8.7 

X 

10~ 5 

1.9 

X 

10" 11 

8.5 

X 

10- 10 

1.6 

X 

10- 6 

5.0 

X 

10~ 5 

7.2 

X 

10 -9 

1.3 

X 

10- 8 

8.1 

X 

T7 8 " 

5.4 

X 

10“ 5 

2.9 

X 

IQ ' 10 

7.0 

X 

10- 9 

1.7 

X 

10- 6 

3.0 

X 

10-' 5 

6.0 

X 

10" 9 

1.9 

X 

10- 8 

5.2 

X 

xo - 8 

3.5 

X 

10-5 

7.1 

X 

10- 10 

6.3 

X 

10~ 9 

1.7 

X 

10- 6 

2.6 

X 

10-5 

3.4 

X 

10- 8 

2.8 

X 

10- 8 

1.7 

X 

10~ 7 

3.3 

X 

10“5 

7.0 

X 

O 

r— 1 
1 

O 

i — 1 

9.1 

X 

10~ 9 

1.6 

X 

10- 6 

8.2 

X 

10- 6 

1.0 

X 

10" 7 

2.2 

X 

10- 8 

1.5 

X 

10~ 7 

3.2 

X 

10- 6 

8.1 

X 

io - : -° 

1.1 

X 

10- 8 

5.2 

X 

10~ 7 

1.2 

X 

10- 6 

1.3 

X 

10” 7 

2.7 

X 

10” 8 

1.0 

X 

10~ 7 

1.1 

X 

10" 6 

5.7 

X 

10~ 9 

5.-4 

X 

10- 8 

6.6 

X 

xo - 8 

3.0 

X 

10“ 5 

2.8 

X 

10- 11 

9.9 

X 

o 

1 — t 
1 

o 

r— 1 

3.4 

X 

xo - 10 

2.9 

X 

10~ 5 

2.8 

X 

10" 11 

1.0 

X 

10~ 9 

4.1 

X 

10 -7 


€ 
















































Table I (Continued) 


ORDER RANK 


30 x 10 2 


30 x 10 5 


30 x 30 


30 x 10 10 


LO x 30 1 


10 x 30 2 


10 x 30 5 


10 x 30 10 


30 x 30 2 


| NORM 1 

NORM 2 

NORM 3 

NORM 4 

2.0 x IO” 8 

2.9 x 10” 9 

1.6 x IO" 9 

2.1 x IO" 8 

1.8 x IO* 5 

1.5 x 10" 10 

2.9 x IO" 9 

6.2 x IO" 7 

X. 4 x IO” 5 

3.0 x 10” 9 

2.0 x IO" 8 

1.8 x IO" 8 

1.5 x IO” 5 

2 . 4 x 10 -30 

2.8 x IO" 9 

4.1 x 10“ 7 

7.5 x IO” 7 

2.8 x IO” 30 

1.4 x IO" 8 

3*5 x IO" 9 

5.4 x IO” 7 

2.2 x IO” 30 

5.8 x IO -9 

5.0 x IO" 9 

7.2 x 10-5 

4 . 9 x 10‘ 3? 

3.5 x IO -30 

3.8 x IO -10 

9.8 x 10“5 

1.1 x IO” 33 

1.1 x IO" 9 

3.4 x IO" 6 

3.5 x lO" 5 

2.1 x IO” 8 

2.4 x IO -8 

1.2 x IO" 8 

5.0 x 10 ” 5 

1.7 x 10" 30 

4.6 x IO -9 

3.8 x IO" 6 

3.4 x 10~5 

4.8 x 10" 8 

4.2 x IO -8 

4. ,3 x IO" 7 

6.0 x lO" 5 

1.3 x 10 -9 

1.1 x IO -8 

5.6 x IO” 6 

6.5 x IO” 7 

2.6 x iO -30 

1.9 x IO" 9 

7.7 x IO” 9 

ij 9.1 x 10" 7 

3.4 x 10" 10 

1.7 x 10“ 8 

7.1 x IO” 7 

. _/i 

1.6 x 10 

1.7 x 10" 13 

7.2 x IO" 10 

3.1 x 10” 1(3 

-=r 

i 

o 
1 — 1 

X 

VO 

• 

rH 

1.8 x IO -13 

5.8 x IO" 10 

4.2 x IO” 6 

• -ii 

2.4 x 10 

1.7 x 10" 8 

1.9 x IO -8 

5.8 x IO” 8 

2.6 x 10 

1.4 x 10 -10 

2.0 x IO -9 

4.3 x IO” 6 


M 



























































Table I (Continued) 


ORDER RANK 


30 x 30 5 


30 x 30 10 


30 x 30 15 


30 x 30 20 


30 x 30 30 


35 x 35 1 


35 x 35 2 


35 x 35 5 


35 x 35 10 


1.2 x 10 
1.2 x 10 

1.0 x 10 

1.1 x 10 


L 

NORM 2 

NORM 3 

NORM 4 

r 11 

9.8 x 10 -9 

2.5 x 10“ 8 

7.2 x 10 

r 1 ' 

3-5 x 10~ 10 

9.6 x 10“ 9 

3.7 x 10 


2.0 x 10" 8 

2.9 x 10“ 8 

2.2 x 10 

r* 1 

7.4 x 10~ ]0 

1 ^ 

• 

I »— • 

X 

o 

I 

oo 

4.8 x 10 

r 8 ' 

8.6x3 0~® 

4.2 x 10“ 8 

3.3 x 10 

)“5 

1.4 x 10“ 

1.1 x 10“ 8 

5.3 x 10 


1.4 x 10“ 7 

3.5 x 10“ 8 

7.7 x 10 


1.2 x 10“ 9 

1.2 x 10" 8 

2.8 x 10 


3.0 x 10 ^ 1.4 x 10 

2.9 x lO" 5 1.2 x 10' 

1.3 x 10“ 6 3.5 x 10' 

1.4 x 10" 6 3.0 x 10' 

3.9 x 10 _l1 1.5 x 10' 

3.8 x 10~ J| 1.5 x 10' 

2.6 x 10" 11 7.9 x 10' 

2.7 x 10" 1 ' 1.5 x 10' 


-11 


-11 


-10 


1.8 x 10' 


2.0 x 10' 


1.5 x 10' 

2.6 x 10' 


1.8 x 10"° 
6.0 x 10“ 10 


3.5 x 10“ 6 
8.1 x 10“ 10 


1.7 x 10 u 

4.7 x 10“ 8 

8.7 x 10 -10 
4.9 x 10“ 10 

1.2 x 10“ 8 

6.2 x 10“ 9 

2.4 x 10“ 8 

9.3 x 10“ 9 


7.5 x 10“ 8 


2.6 x 10 
7.9 x 10“ 6 

1.7 x 10~ 7 
5.2 x 10“ 6 

2.1 x 10 

6.7 x 10 


2.7 x 10“ 7 
6.0 x 10' 6 


t 


• • 








Table I (Continued) 


ORDER RANK 


35 x 35 15 


35 x 35 20 


35 x 35 30 


35 x 35 35 


x 15 



*10 x 15 ? 


*10 x 15 5 


*10 X 15 10 


*10 x 15 15 


NORM 

1 1 

NORM 

1 2 

NORM 

1 3 

NORTd 

[ 4 

1.5 

X 

io“ n 

6.1 

X 

10 - 8 ' 

* 1.5 

X 

10- 8 

4.7 

X 

10~ 7 

1.6 

X 

10“ 11 

7.3 

X 

10- 10 

1.0 

X 

10- 8 

5.1 

X 

io - 6 

2.4 

X 

10“ 5 

1.2 

X 

10" 7 

3.3 

X 

10- 8 

3.1 

X 

10“ 7 

1 . ] 

X 

10~ 5 

8.8 

X 

10- 10 

1.2 

X 

10“ 8 

4.0 

X 

io - 6 

1.8 

X 

10“ 5 

2.6 

X 

10“ 7 

2.6 

X 

10“ 8 

2.7 

X 

10“ 7 

1.3 

X 

10“ 5 

1.6 

X 

10“ 9 

2.0 

X 

10- 8 

1.5 

X 

io - 6 

1.8 

X 

JO ' 6 

3 . 4 

X 

10“ 9 

2.3 

X 

10- 8 

1.5 

X 

10“ 7 

1.8 

X 

10“ 5 

* 1.0 

X 

10" 9 

6.5 

X 

10- 8 

7.0 

X 

io - 8 

7.7 

X 

10 “ 5 

2.5 

X 

IO " 11 

9.0 

X 

IO -10 

4.0 

X 

io - 10 

7.8 

X 

10- 5 

2.5 

:c 

10- 11 

6.5 

X 

io - 10 

1.2 

X 

io - 6 

7.0 

X 

10~ 5 

6.5 

X 

10“ 9 

1.8 

X 

ID " 8 

6.4 

X 

] 0“ 8 

7.1 

X 

10' 5 

1.6 

X 

10- 10 

3.1 

X 

10“ 9 

5.1 

X 

io - 7 

5.8 

X 

10“ 5 

3.9 

X 

10“ 9 

2.3 

X 

io - 8 

3.6 

X 

io - 8 

5.6 

X 

10“ 5 

2 . *1 

X 

10- 10 

4.6 

X 

10“ 9 

4.3 

X 

10“ 7 

2.1 

X 

10” 5 

5.3 

•X 

10-9 

2.0 

X 

io - 8 

1.9 

X 

io - 8 

2.0 

X 

10“ 5 

2.2 

X 

10- 10 

4.8 

X 

IO ” 9 

5.1 

X 

io - 7 

8.9 

X 

10- 7 

2.7 

X 

10- 10 

1.4 

X 

io - 8 

3.1 

X 

10“ 9 

9-0 

X 

10“ 7 

2.9 

X 

♦— » 
o 

1 

1 —* 

o 

6.7 

X 

10“ 9 

6.4 

X 

IO " 9 

























































Table I (Continued) 


OHDEK HANK 


40 x 40 



NORM ] 


4.4 x 10“ 


NORM 2 


1.2 x 


4.3 x 10 1.2 x 


NORM 3 


6.0 x 
4.8 x 


,-10 


,-10 


NORM 4 



40 x 40 2 


4.1 x 10' 


1.2 x 


4.] x 10 9.2 x 


1.6 x 


5.7 x 


10~ 9 


1.6 x 10' 
9.4 x 10* 


40 x 40 5 


2.0 x 10' 


3.2 x 10' 


3.1 x 
1.6 x 


,-10 


7.4 x 


1.1 x 


1.7 x 10' 
7.0 x 10’ 


40 x 40 10 


1.4 x 10' 


3.0 x 


1.9 x 10 6.5 x 


9.1 x 
1.4 x 


8.5 x 10’ 


40 x 40 15 


1.2 x 10 


2.1 x 10' 


2.6 x 

9.7 x 


,-10 


8.1 x 
1.3 x 


10- 8 


7.2 x 10' 


40 x 40 20 


1.3 x 10" 6.7 x 


1.4 x 10' 


1.0 x 


4.3 x 


1.1 x 


4.7 x 3 0’ 


40 x 40 30 


3.7 x 10*° 6.2 x 


3.5 x 10' 


1.2 x 


3.0 x 


1.2 x 


5.3 x 10' 
4.5 x 10' 


40 x 40 35 


40 x 40 39 


1.1 x 10 


5.5 -x 10 


1.0 x 10~ 5 1.9 x 10' 


5.4 x 10"° 1.1 x 10' 


7.2 x 10' 


2.1 x 10 


3.1 x 
2.7 x 


3.8 x 


1.0 x 10~ u 
1.3 x 10~ 6 


2 .If x 10~ 7 
3.8 x 10" 7 
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Table I (Concluded) 


ORDER 

RANK 

NORM 1 

NORM 2 

NORM 3 

NORM 4 

40 x <10 

40 

2.6 x 10" 6 
2.5 x 10" 6 

7.1 x 10" 9 

8.2 x 10" 9 

3.7 x 10" 8 
9.3 x 10" 8 

2.6 x 10" 7 
1.3 x 10" 7 

45 x 

30 

7.2 x 10 -5 
6.9 x 10-5 

6.7 x 10" 8 
7.9 x 10 _1 ° 

3.0 x 10" 8 

1.1 x 10“ 8 

2.8 x 10" 7 
3.7 x 10" 6 

<15 x <10 

35 

2.2 x 10”5 
1.8 x 10“ 5 

1.2 'x 10~ 7 
9 . <* x 10" 10 

2.8 x 10“ 8 
1.4 x 10" 8 

1.'5 x 10- 7 
3.1 x 10“ 6 

<15 x <10 

39 

2.5 x 10“ 5 
5.7 x 10" 6 

<1.8 x 10" 7 
9.0 x 10“ 10 

1.9 x 10 -8 
1.7 x 10" 8 

1.3 x“IO' 7 
1.3 x IQ' 7 
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Table II 


Computer Times Required to Calculate the Generalized 
Inverse of a 10 x 10 Nonsingular Matrix 


PROGRAM 

MEAN TIME (MILLI-SECONDS) 

G1NV2 

3*1.9 

SEQINV 

50.3 

JPLUS 

390.5 

APLUS 

1117.2 


Table III 

Computer Times Required by GINV2 and APLUS 
for Selected Randomly Generated Matrices. 


MATRIX 

MEAN TIME (MILLI-SECONDS) | 

• 

TYPE 

ORDER 

GINV2 

APLUS 1 

NONSINGULAR 

10 x 10 

3H.9 

417.2 

NONSINGULAR 

20 x 20 

245.8 

3391 * 

•NONSINGULAR 

30 x 30 

784.1 

10981 ** 

NONSINGULAR 

40 x 40 

1825.0 

27383 *** 


(Note: The mean is, in general, calculated for a block of 

100 matrices) 

^Calculated for a block of only 50 matrices 

**Calculatea for a block of only 25 matrices 

***Calculat ed for a block of only 15 matrices 
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Table IV 

Computer Times Required by GINV2 and APLUS 
for Selected Randomly Generated Matrices. 


MATRIX 

ORDER 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

X 

10 

10 

• 

X 

10 

20 

X 

20 

20 

X 

20 

20 

X 

20 

20 

X 

20 

20 

X 

20 

20 

X 

20 

20 

X 

20 


RANK 



MEAN TIME ( MI LLI -SECONDS) 

GINV2 

APLUS 

43 

29 

43 

377 

42 

395 

42 

404 

41 

4l4 

40 

436 

39 

442 

38 

46l 

37 

477 ' 

35 

4 16 

301 

266 

300 

3243 

300 

3204 

299 

3354 

294 

3254 

293 

3390 

291 

3409 













Table IV (Continued) 


MATRIX 

TYPE 

MEAN TIME (MILLI-SECOHDS) 

; ORDER 

RANK 

GINV2 

APLUS 

! 20 x 20 

8 

287 

31468 

• 

20 x 20 

9 

283 

3 H 61 

20 x 20 

10 

282 

3556 

20 x 20 

11 

282 

3533 

20 x 20 

12 

278 

3612 

20 x 20 

13 

277 

3730 

•20 x 20 

lH 

275 

3776 

20 x 20 

/ •* 

15 

269 

379'i 

20 x 20 

16 

26H 

3710 

20 x 20 

17 

259 

39 H 3 

20 x 20 

18 

25 ^ 

3857 

20 x 20 

19 

2 kg 

3993 

20 x 20 

20 

2k3 

3354 
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APPENDIX A 


THEORETICAL BACKGROUND 

The following is a brief ..description of the mathematics 
used in each of the five programs examined in this study. 

<1) PEN 2 is based on a method, devised by R, Penrose 
♦ 

(ref. 2), which computes the generalized inverse 
of a matrix A using the formula 

(a) A + = DA T 

where D is any matrix satisfying 

(b) A T A = D(A T A) 2 

/ \ 

| 

(Note that multiplication of (a) on the right 
by A + (A + )V yields (b).) 

Define a sequence of matrices C^, k = 1,2,... 

by ' 


Ci = i 


C 2 = I • Y trace (CxB) - CiB 


C, X1 = I • r trace (C.B) - C.B 
j+i j j J 


where B = A A. 


If r is the rank of A, then C ,.B = 0 and 

9 r+ 1 


trace C r B ^ 0. 


• A-l 


#* 


I » 

* 

* » 

Then D can be calculated by the formula 


(2) SEQ1NV is based on a sequential method for com- 
puting the generalized inverse (ref. 9)* Let A n-1 

be the matrix containing the xirst (n~l) rows of 

* th 

the matrix A and a n be the n row of A. The 

generalized inverse of A is then calculated 

, sequentially by the formula 



p a A 
^n n n-i 



with p n = 


I* # 

^a n (I - A n -i A n-i^ a n ^ a n A n-i A n-i 

fl + a (A + ) T a T V 1 A 4 " (A + ) T a T if a 
y n n - 1 n J n-i K n-i' n n 


aA M , A . 
n n-i n-i 


j. .j. m ' 

A n A n anc * A n^ A n^ are c 9 m P u ^ eci sequentially as 
follows: 


A + A 
♦ n A n 



- p M a n A M , 
^n n n-i 



4- 

'-A A -paA ,A ,+pa 
n-in~i p n n n-i n-i p n n 


If a n = a n A n-i A n-i’ then A n A n “ A n-i A n-i 


A-2 


then k+k 
n n 


A n-i A n~i + 'P n Pn • 


If a M / a A A 
n n n-i n- 



p a A 
n 


+ 

n-i 



” A n-i( A n~i) T " p n a n A n-i( A n-i) T 
+ Pn( a nCi'C/ a I + X ) p n 



(3) The method JPLUS uses the property that if A is an 

\ 

n x n symmetric matrix, there exists an n x »n 
matrix P such that . • 

T 

PAP 1 = D 


where D is the diagonal matrix whose diagonal 
elements are the eigenvalues of A. It can be 
shown that 

« 

+ T + 

A = P A D P . 


A-3 



If A is an n x m matrix, then let B = A A, Then, 
by the above, there is an m x m matrix P such that 

T 

D = PB? 1 

* , 

where D is again diagonal with the eigenvalues of 
B as its diagonal elements. 

r 

Then 

& _+ Tt-, 'f'rs 

B = PDP . 


Using the fact that 


A + = '(A T A) + A , 


it follows that 


A + = B + A. 


(*l) APLUS in an iterative method (ref. 10) based on 

» 

the following formula 


X n = X n-x < 21 " AX n-i> 


where A is the matrix whose generalized inverse is 
to be computed. After, initially setting 

m .... 

a 1 ■ , 

X - — the iteration process is continued 

0 1 1 aa t 1 l 

until a near-zero test indicates that X^ can be 
assumed to be A + for some m, (Note that here, 


• / n m / m \ \ 

1 1 aa t 1 1 = ( E E ( AA )t; ) 

\i=l j=l\ /- J / 


Vz 


.) 


A -4 




* 

(5) GINV2 is based on an extension of the orthogonal- 
ization process for computing the Inverse of non- 
singular matrices (ref. 11). The problem of 
computing the generalized inverse of an n x m 
matrix A can be reduced to the problem of com- 
puting the inverse of an n x m matrix (R*S), 
partitioned so that R is the matrix of all 
linearly independent columns (say there are k of 
them) of A and S is the matrix consisting of the 
remaining (m-k) dependent columns. (Note that 
(Rjs) can be obtained from A by a finite number 
of permutations of columns of A.) 


A Gram-Schmitt orthogonalization process is per- 
formed on (R*S) and the same operations performed 
simultaneously on the n x n identity matrix. The 
result of the Gram-Schmitt on (R;S) is a matrix 
of the form (Q*0) where Q is an n x k matrix and 
0 is the n x (m-k) zero matrix. The identity 
matrix will become of the form 




where Z is k x k 


0 is the (n-k) x k zero matrix, ~U is k x (n-k) 
and I . is the (n-k) identity matrix. 

li"* - K * * 

Performing a Gram-Schmitt orthogonalization on 

\ 

j yields a matrix of the 
A-5 



i 



form 



where P Is (n-k) x (n-k) and 


«UP is k x n-k. 


After these orthogonalizations are completed, all 
of the matrices necessary for the computation of 

4 . 

A have been calculated. 


It can be shown that A 


is given by 




q t 

l(UP) T ZQ T 



- (UP)(UP) T ZQ T 

p(up) t zq t 


APPENDIX B 


ANALYSIS OP VARIANCE 


This appendix describes the analysis of variance used 
to test the hypothesis that the accuracies of GINV2 and APLUS 
are equal, The observations used in the analysis were the 
absolute values of the logarithms of the norms, which were 
defined on page 6 of the report. 

CASE I. Full rank. 

A three-way full factorial was assumed as the model, i.e. 


ijk 


M + a ± + 0^ + + ((*3).^ + (“^ik + (^) jk + P 


ijk 


where : 




i = 1,2 

j - 1 , 2 , . . . , 15 
k = 1,2, 3, 4 

a x : effect due to APLUS 
a 2 : effect due to GINV2 

i » 

■ Bj : effect due to order (orders considered were— 
(2 x 2), (2 x 3)/ (3 x ' 2 ) , '(4 x 5), (10 x 8), 
: (8 X 10), (10 X 10), (20 X 20), (30 X 10), 

(10 x 30), (30 x 30), (35 x 35), (40 x 15), 
(40 x 40), (45 x 40)) 


B-l 




* * 



effect 

due 

to norms 






(op) ij : 

effect 

of 

interaction 

due 

to 


and 

6 j 

( aTf) jk : 

* effect 

of 

interaction 

due 

o 0 

“i 

and 

Y k 

^ 6Y ^jk : 

effect 

of 

interaction 

due 

to 


and 

Y k 

p ijk : 

effect 

due 

to random < 

error 





The null hypothesis is as follows: 


H q : accuracy of GINV2 « accuracy of APLUS 
The results of the analysis are tabulated below: 


Analysis of Variance' 

Source of 


Variation 

D.P. 

• SS ’ . 

• * MS 

P 

Method 

1 

0.01391 

0.01397 

0.1405 

Order 

14 

19.3749 

1.3837 

13.91 

Norm 

3 

94.9728 

31.6576 

318.89 

Method x Order 

14 

7.0123 

0.5009 


Method x Norm 

3 

1.-4788 

0.4930 


Order x Norm 

42 

12.2357 

0.2912 


Error 

42 

4.1762 

0.0994 



Since the P ratio for the methods is .1405 it is not 
possible to reject H q based on the data so far collected. 


B-2 


CASE II, Nonfull rank. * 

The model assumed In this case was a four-way 
partially-nested factorial, in which rank is nested in any 
order. The results of this case' are tabulated below. 

H : accuracy of GINV2 = accuracy of APLUS 

o 

t 

Analysis of Variance 


Source 

of Variation 

D.’P. 

ss 

MS 

P 

Method 


1 

2.596 

2.596 

60 

Order 


9 

22.525 

2.503 


Hank 


1 

0.633 

* 

0.033 


Norm 

* 

3 

324.595 

108.198 


• 

• 

• 

Error 

Interactions 

• 

• 

• 

• 

• 

• 

I.1496 

• 

• 

• 

0.0425 



Since the P ratio for the method is 6 0, the hypothesis 
that the two methods are equal must be rejected in favor of 
the hypothesis that the accuracy of GINV2 is better than 


that of APLUS. 


